A model for the epitaxial growth of graphene on 6H-SiC 
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We introduce a kinetic model for the growth of epitaxial graphene on 6H-SiC. The model applies 
to vicinal surfaces composed of half-unit-cell height steps where experiment shows that step flow 
sublimation of SiC promotes the formation and growth of graphene strips parallel to the step edges. 
The model parameters are effective energy barriers for the nucleation and subsequent propagation 
of graphene at the step edges. Using both rate equations and kinetic Monte Carlo simulations, 
two distinct growth regimes emerge from a study of the layer coverage and distribution of top-layer 
graphene strip widths as a function of total coverage, vicinal angle, and the model parameters. 
One regime is dominated by the coalescence of strips. The other regime is dominated by a novel 
"climb-over" process which facilitates the propagation of graphene from one terrace to the next. 
Comparing our results to scanning microscopy studies will provide the first quantitative insights 
into the kinetics of growth for this unique epitaxial system. 

PACS numbers: 81.15.Aa, 68.55.-a, 68.35.-p 
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I. INTRODUCTION 

Of the various ways known to produce graphene, it is 
increasingly likely that only graphene grown epitaxially 
on silicon carbide will play an important role in post- 
CMOS microelectronics^. This system is unique among 
epitaxial growth systems because there is no deposition 
flux. Instead, silicon atoms sublime from SiC(OOOl) and 
SiC (0001) at high temperature and the carbon atoms 
left behind recrystallize into graphene. Microscopy stud- 
ies of vicinal 6H-SiC surfaces show that graphene nu- 
cleates at step edges and that the subsequent morphol- 
ogy depends strongly on growth conditions, vicinality, 
and whether the step heights are equal to one, two, or 
three Si-C bilayers 3 - - — . For example, growth from single- 
bilayer steps produces a complex, finger-like graphene 
morphology 5,8 while step-flow growth from triple-bilayer 
(half-unit-cell) height steps produces long, straight strips 
parallel to the steps whose widths increase as growth 
proceeds 3,8 . The difference comes from the fact that 
three bilayers of SiC must desorb to liberate a sufficient 
number of carbons atoms to cover the sublimed area with 
one layer of graphene. 

The calculations reported in this article aim to (i) pro- 
vide experimenters with a simple and convenient way to 
characterize the changes they see in surface morphology 
when growth conditions change; (ii) identify a statis- 
tical measure of sub-monolayer growth which identifies 
whether graphene step-flow growth is limited by nucle- 
ation at steps or by propagation on terraces; and (iii) 
provide physical insight into the competition between 
graphene strip coalescence and a new kinetic process 
(unique to this system) which we call "climb-over" . Our 
principal theoretical tool is a phenomenological kinetic 
Monte Carlo (KMC) model of the sort used to study the 
growth kinetics of III-V semiconductors^. When cou- 
pled closely with experiment, this approach produced a 
decade of valuable insights before simulations based on 
total energy calculations of energy barriers for III-V sys- 



tems became possible^. For our problem, first-principles 
KMC is impossible because the structure of the "buffer 
layer" at the graphene/SiC interface is controversial and 
the structure of steps on this buffer layer is unknown 11 . 
On the other hand, the fact that graphene grows in strips 
from triple bilayer steps means that a one-dimensional 
model is a good first approximation and the distribution 
of these strip widths will play a prominent role in what 
follows. A mean- field rate-equation analysis in the Ap- 
pendix provides further insight into our proposed model 
and the KMC results. 



II. KINETIC MONTE CARLO MODEL 

Fig. [1] shows the various processes we consider for a 
vicinal surface of 6H-SiC composed exclusively of half- 
unit-cell height steps. Each process involves the replace- 
ment of a unit area of SiC triple bilayer by a unit area 
of graphene. Consistent with the coarse-grained nature 
of the our model, we do not concern ourselves with 
atomistic details and simply assume that all exposed 
SiC terraces spontaneously reconstruct to the carbon- 
rich "buffer layer" known to form on both 6H-SiC(0001) 
and 6H-SiC(0001) u . Fig. [Ha)-(b) shows the nucleation 
of graphene at a SiC step with no graphene nearest 
neighbors. This occurs in our model at a rate r nuc = 



v exp(-E nuc /kT), where v 



10 1 



is an attempt 



frequency and T is the substrate temperature. Two 
points are worth noting. First, E nuc is an effective en- 
ergy parameter which accounts for the combined effects 
of Si atoms sublimation, C atoms re-crystallization, and 
graphene growth along the step edge. Second, a variation 
of our model could allow additional SiC to sublime be- 
fore a stable graphene nucleus forms. This influences the 
predicted distribution of strip widths and, like the corre- 
sponding problem of critical island sizes in conventional 
epitaxial growth, comparison with experiment provides 
microscopic information that is nearly impossible to learn 
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FIG. 1. Kinetic processes allowed in the KMC simulation. 
ThMK6DTY-T7S9SCe steps marked A (green), B (red), C 
(blue), and D (purple) play a role in the rate theory reported 
in the Appendix. 



any other way*2. 

After nucleation, graphene growth continues by dis- 
solution of the adjacent SiC step at a rate r prop = 
1/q exp(— E prop /kT) [Fig. [T|b)-(c)]. Propagation occurs 
only at SiC steps that are bounded by a graphene strip. 
Two fates are possible for such a strip. One is that the 
propagating strip runs into another SiC step and creates 
a step bunch of two triple bilayer steps [Fig. Hfd))]. If 
this happens, the strip can "climb over" the upper ter- 
race at the rate r prop [Fig. Hfe)]. Another possibility is 
that the propagating strip meets another strip on the up- 
per terrace [Fig.QTf)]. In this case, our KMC simulation 
coalesces the two strips at the rate r prop [Fig. QJg)]. Nu- 
cleation of a covered graphene layer at a covered SiC step 
occurs at the rate r nuc [Fig. njg)-(h)]. Propagation of a 
covered graphene layer occurs at the rate r prop or (for 
some of the simulations reported below) at the slower 
rate r prop [Fig. HJi)]. The later growth continues in the 
same way as the first graphene layer. 

We use a standard KMC algorithm 13 to simulate 
growth on vicinal SiC surfaces composed of (typically) 
5000 steps with periodic boundary conditions. At least 
100 independent runs were averaged to obtain sta- 
tistically significant results. The vicinal angle (f> — 
tan _1 (3/iy), where W is the terrace width. We be- 
gin our discussion with 0^, the graphene coverage of 
layer i, as a function of the total graphene coverage 
= These quantities are accessible to spatial- 

averaging experimental probes and our model energy pa- 
rameters should provide a simple and convenient way for 
experimenters to characterize variations in observed mor- 
phology with growth conditions. Later, we will turn to 
the distribution of graphene strip widths as a quantity 
which scanning microscopy can exploit to learn the rel- 
ative importance of competing surface kinetic processes 
during growth. 



AE/kT 

FIG. 2. Layer 1 coverage as a function of the energy bar- 
rier difference AE and vicinality (p. Solid and dashed lines 
correspond to (f> = 0.9° and <j) — 3.4°, respectively. 



III. LAYER COVERAGE AND GROWTH TIME 

Fig. [2] shows simulation results for ©i as a function of 
AE = E rmc - E pmp with r prop = r prop for four different 
values of total coverage and two choices for the vicinal 
angle <f>. The rather counter- intuitive behavior that 0i 
decreases as AE increases for fixed can be understood 
as follows. When AE is large, propagation of existing 
graphene strips is relatively more likely than the nucle- 
ation of new graphene strips, and fewer graphene strips 
can form. Many strips undergo the "climb over" pro- 
cess [Fig.[ljd)-(e)] when the width of the graphene strips 
passes the terrace width, thereby creating many nucle- 
ation sites for second layer growth. The net result is that 
nucleation of second layer graphene begins earlier. Be- 
cause these strips grow for a longer time, the total second 
layer coverage is larger. Fig. [2] also implies that better 
surface homogeneity can be achieved by increasing the 
substrate temperature and decreasing the substrate mis- 
cut angle. This conclusion is consistent with the observa- 
tions reported in Ref. [IJ. When AE is further increased, 
the number of nucleation events for both the first layer 
and second layer are greatly reduced, and eventually the 
competition between the two layers is balanced. This 
effect produces the lower plateaus in Fig. [2] A rate equa- 
tion analysis described in the Appendix provides another 
way to understand this behavior of our model. 

In principle, experimental data for and 0i can be 
compared with the curves in Fig. [2] (outside the plateau 
regime) to extract a value for AE. However, because it 
is surely harder for Si atoms to escape from SiC when 
they are covered by a graphene layer than when they 
are not, we introduce a second layer propagation barrier 
-^prop > -^prop- The corresponding rate for second layer 
propagation is r prop = uq exp(— E prop /kT) and we define 
AE' = E' prop - E plop . We retain the equality of the first 
and second layer nucleation barriers for simplicity^. We 
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FIG. 3. Layer 1 coverage as a function of AE for differ- 
ent values of AE' . AE' /kT = 0,0.6, 1.3, 1.9 applies to the 
solid curves, dashed curves, dashed-dotted curves, and dot- 
ted curves, respectively. The vicinal angle <j) = 0.9°. 
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FIG. 4. KMC time tk as a function of the energy barrier 
difference AE and AE' at a fixed total coverage. 



also forbid the growth of layer 3. Fig.f3]shows that as AE' 
increases, the first layer coverage increases substantially, 
as might be expected. We now have a three parameter 
problem and experimental data for 8i at two values of 
total coverage O can be used to extract values for AE 
and AE' from Fig. [3] 

It remains to deduce values of -E^uo E piop and E' 
individually. This can be done using the experimental 
growth time for a given total coverage because tk = 
tEfpmp relates the dimensionless KMC simulation time 
to the experimental growth time tE- Fig. |4] shows tk 
as a function of AE for different values of AE' . This 
graph (or a similar one obtained for a different choice of 
and 4>) permits -Eprop to be extracted from the values 
of AE and AE 1 determined earlier from the layer cov- 
erage curves. The two other energy parameters follow 
immediately. 




FIG. 5. (a) LEEM image of graphene grown on vicinal 6H- 
SiC(0001) fromRef.0. Re gions covered by one, two, and three 
layers of graphene are shown as light, moderate, and dark 
gray, respectively. The latter two occur at SiC step edges. 
(b)-(d), KMC simulation images of monolayer graphene strips 
with AE/kT — 0, 5.8, and 11.6, respectively. The total cover- 
age Q = 0.25. Light grey lines and the right edges of graphene 
stripes are SiC steps. The vicinal angle <f> = 0.9°. 



IV. STRIP WIDTH DISTRIBUTION 

We turn now to the distribution of first-layer graphene 
strip widths. This statistical quantity probes more 
deeply into the competition between nucleation and prop- 
agation and between coalescence and climb-over. It also 
provides another way to extract AE and to understand 
the cross-over from the \ow-AE plateau to the high-Ai? 
plateau in Fig. [2j Compared to the coverage curves, this 
distribution is much more sensitive to AE and much less 
sensitive to AE'. For that reason, we set the latter equal 
to zero in what follows. 

Fig. [SJa) shows a LEEM imaged where the terraces 
are mostly covered by a single monolayer of graphene 
(light gray). Very near the step edges, strips composed 
of two (moderate gray) and three (dark gray) layers of 
graphene are apparent. Fig. EJb) shows a KMC simu- 
lated morphology (with AE = 0) which looks quite simi- 
lar to Fig.[SJa). The graphene strip morphology changes 
significantly as AE increases in Fig. [SJc) and (d): the 
number of graphene strips decreases and many of them 
cover many SiC steps. Note also the change in scale 
from Fig. To quantify this morphological change, 

Fig. [5] plots p(s), the normalized distribution of strips 
with width s, for different choices of AE and O. The 
terrace width here is W = 200. When AE/kT = [Fig- 
ure [BJ a)], the distribution is Poisson because graphene 
nucleates at almost all the SiC steps simultaneously. The 
mean strip width is WO in the interval [0, W]. However, 
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persists when the value of AE increases. However, a 
larger value of AE implies that some graphene strips nu- 
cleate earlier than others. This leads to a shift to the 
right in the peak position seen in Fig. [Bfa)-(c) for the 
same coverage. The increasing time delay between con- 
secutive nucleation events similarly produces a distinct 
broadening of the distribution curves. Eventually, for 
large enough AE, the distribution curves become uni- 
form [Fig. [6fd)] in the scale we consider. This occurs 
when the graphene strips propagate so rapidly (relative 
to nucleation) that the step edges are no longer distin- 
guishable. 

Finally, we return to the vicinal angle dependence of 
the coverage curves plotted in Fig. [2] The p(s) results 
above imply that the transition between the two hori- 
zontal plateaus in these graphs as AE increases reflects 
a transition from a Poisson distribution to a uniform dis- 
tribution of graphene strip widths. In the Poisson regime, 
the terrace width only affects the coverage distribution 
at late times when the maximum graphene strip width 
passes W. Therefore, a change in the vicinal angle 4> 
only changes the coverage distribution for large 0. This 
may be contrasted with the uniform regime, where the 
graphene strips grow so rapidly that they are not hin- 
dered by the SiC step edges. In this case, the coverage 
distribution does not depend on <j> at all. Nevertheless, 
as we see from Fig. [5J for fixed 0, 0i tend to be larger 
for smaller vicinal angle <j>. This is so because the stan- 
dard deviation divided by the terrace width for a Poisson 
distribution is y/WQ/W ~ 1/y/W, which implies that a 
smaller vicinal angle leads to a relatively narrower distri- 
bution of strip widths. In the limit when all the graphene 
strips are about the same width, coalescence events occur 
only very near = 1 and there is essentially no second 
layer growth. This supports our previous statement that 
better surface uniformity can be achieved by using a more 
singular surface. 



CONCLUSION 



FIG. 6. Graphene strip width distribution p(s) for different 
AE and total coverages with <j> = 0.9°. Different color lines 
correspond to = 0.1 (red), 0.3 (magenta), 0.5 (green) and 
1.0 (blue), respectively. The terrace width is W — 200. 



when the coverage > 0.8 (blue line in Fig. Eta)), the 
leading edge of the Poisson distribution begins to cross 
the terrace width W , a few strips disappear by the "coa- 
lescence" mechanism [Fig.[TJf)-(g)], and a few strips with 
widths close to 2W form. As a result, the strip width dis- 
tribution is abruptly cut-off at the terrace width and the 
Poisson distribution repeats (with a much decreased peak 
amplitude) in the width interval [W, 2W]. The distri- 
bution moves farther across the terrace width boundary 
when increases further. 

The general behavior of p(s) with increasing coverage 



In summary, we have developed a one-dimensional ki- 
netic Monte Carlo model to study the epitaxial growth 
of graphene by the step flow sublimation of vicinal sur- 
faces of 6H-SiC(0001). The layer coverages and the dis- 
tribution of graphene strip widths were found to depend 
more or less strongly on the relative sizes of the effec- 
tive energy barriers for graphene nucleation, first layer 
propagation, and second layer propagation. The cross 
over of the distribution from Poisson to uniform as the 
nucleation barrier increases clearly shows that there are 
two distinctive growth regimes, one dominated by "coa- 
lescence" processes and one dominated by "climb over" 
processes. The "climb over" processes have the effect of 
increasing the graphene surface inhomogencity. It will be 
interesting to compare these simulation results with ex- 
perimental measurements to see how the effective energy 
barriers depend on growth parameters like the partial 



pressure of silicon in the growth chamber. Future sim- 
ulations studies will examine the "kinetic roughening" 
of this new model of epitaxial growth^, and the ability 
of the model to rationalize the non-uniform layer thick- 
nesses observed when graphene grows on spontaneously 
facetted SiC substrates^. 
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Total Coverage, 

FIG. 7. The second layer coverage ©2 as a function of the 
total coverage with AE' — 0. The solid lines are KMC 
simulations with (bottom to top) AE/kT = 0,3.9,5.8 and 
11.6. Dashed Lines are the rate equation results. 



VII. APPENDIX: RATE EQUATION ANALYSIS 

This Appendix presents a mean-field rate-equation 
analysis to provide further understanding of our graphene 
growth model. The color coding in Fig. [T] identifies four 
basic types of steps: (A) a bare SiC step; (B) a step con- 
nected to a layer 1 graphene segment; (C) a step that 
is carpeted by a continuous layer of graphene; and (D) a 
step that is connected to a layer 2 graphene segment. We 
let riA, tib, nc, and ri£> be the number of steps of each 
type, so n = ua + + nc + njj is the total number of 
steps and L — nW is the system size. Then, if pans is 
the number of B-steps with an A-step immediately above 
[Figure |TJd)] and Pfn B is the number of B-steps with 
a graphene segment immediately above [Figure [TJf)], an 
approximate description of the epitaxial graphene growth 
processes is 



driA 
~dB 
dns 
~d& 

dn c 

de 

dnp 
~d® 



-fin a ~ r 2 p d n B 



r\n a - r 2 pfn B 



-r\n c + r 2 (pd + Pf)n B 



(1) 
(2) 
(3) 
(4) 



r P To P L/r tot , and r tot 



where n = r DUC L/r tot , r 2 

nc + r' prop n D . 
Eq. (fTJ) says that A-steps (green) are lost by first-layer 
nucleation events and by climb-over events. Eq. ([2]) says 
that B-steps (red) are created by nucleation events and 
lost by coalescence events. Eq. ((3|) says that C-steps 
(blue) are lost by second-layer nucleation events and cre- 
ated by both climb-over and coalescence events. Eq. (|4]) 
says that D-steps (purple) are created by second-layer 
nucleation events. We note that a climb-over event does 
not change the number of B-steps. 



We consider two limits where pd and pf can be esti- 
mated. The first limit is r nuc -C r prop where first-layer 
nucleation events are rare. Climb-over is frequent and 
coalescence infrequent. These conditions imply, in turn, 
that pf <C 1 and pd ~ 1/W. The second of these is 
true because, when the coverage is fixed and the prop- 
agation rate is very fast, the length (modulo W) of the 
graphene segment connected to a B-step takes every value 
between one and the terrace length W. Conversely when 
r nuc = r prop , nearly every step produces a nucleation 
event and climb-over is rare. This implies that pd <C 1 
and pf is the probability that the length (modulo W) of 
the graphene segment connected to a B-step is W — 1 
as determined from a Poisson distribution with average 
value WO. 

We have solved Eqs. numerically (assuming 

AE' = for simplicity) in the two limits discussed above 
using the initial conditions 



n A 



n B = n c = n D = 0. 



Using this numerical data, we calculate 



dOi dQi/dt r nuc nA 
dO^ ~ d0 2 /dt ~ 



-1, 



(5) 



(6) 



to 



and use = 81 + 202 to equate the right side of 
dO/d& 2 - 2. 

Fig. [7] compares O2 versus as determined from the 
KMC simulation (solid curves) with the correspondingly 
rate equation results (dashed curves). The agreement is 
quite good when AE/kT = (purple curve). This is the 
no "climb-over" regime where a = Pf /{pf +Pd) = 1. The 
agreement is similarly good when AE/kT is large (red 
curve) if we account for coalescence in the rate equa- 
tions with the choice a = 0.01. The no-coalescence 
curve (a — 0) falls below the a ^ curve because, in 
the rate equations, the presence of coalescence reduces 
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the life time for all first layer propagating graphene seg- 
ments, which reduce the number of competitors to sec- 
ond layer propagation. Because AE/kT is large, there 
are not many graphene segments in the system to begin 



with. Removing some first layer segments by coalescence 
promotes second layer propagation and thus results in a 
larger second layer coverage. 
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